* This version does all regressions with fuel consumption and local pollution from AP2
* but for tracks not in main sample

* location of datafiles
local trackdatafile = "${base_dir}/data/final/CA_tracks_other.csv" 

local analysis_dta = "${base_dir}/data/final/CA_tracks_other.dta"
local reassemble = 1

* location of fuel price and cpi data
local price_file "${base_dir}/data_not_to_share/SPfuelprices/SPfuelprices.csv"
local cpi_file "${base_dir}/data/CPI/CPIAUCSL"
local unem_file "${base_dir}/data/unemployment/UNRATE"

* Assembling *******************************************************************
if `reassemble'==1 {

	* make WFR data
	*do "`base_dir'Stata/CA_tracks/load_vessel_char.do" `wfr_folder' `wfr_date' `wfrdta'
	*do "`base_dir'Stata/CA_tracks/load_vessel_list.do" `vesselList' `vesselList_dta'
	
	do "${base_dir}/analysis_code/CAtracks_assemble" ///
		"`trackdatafile'" "`price_file'" "`cpi_file'" "`unem_file'"
	save "`analysis_dta'" , replace
}
else {
	use "`analysis_dta'" , replace
}

* INITIAL CLEANING *************************************************************
drop if date1<mdy(1,1,2009) 			// getting 2008 dates for some reason!
										// likely tracks on 1/1/2009 and time is extrapolated in some way
									
replace length=. if length>450 // clearing really long boats (probably data error)

drop if inlist(group_agg,"Passenger","Miscel")
drop if inlist(group,"Tug","Dredge")

* dropping PW Sound and Cook Inlet routes
drop if otherAK == 1 

* dropping HI to HI routes
drop if HItoHI == 1 

* if track crosses land too much then drop
drop if km_onland>5

* drop if fuel consumption equals zero
drop if f_cons==0 | f_power==0

* update route_id, to be missing if either port is missing
replace route_id=. if port1==. | port2==.

********************************************************************************
* recovering eca dates
levelsof eca_ind if eca_ind>0 , local(eca_inds)
foreach i of local eca_inds {
	qui sum date1 if eca_ind==`i'
	scalar eca`i'd = r(min)
	scalar eca`i'qy = yq( year(r(min)),quarter(r(min)) ) 
	scalar eca`i'my = ym( year(r(min)),month(r(min)) ) 
}


* Other Cleaning ***************************************************************
* creates main_sample variable for use later
do "${base_dir}/analysis_code/cleaning.do"
	
* predict FCC for vessels without characteristics
do "${base_dir}/analysis_code/predict_FCC.do"	


label var dist "Dist, (km)"
label var f_cons "Fuel, (t)"
label var faux "Aux Fuel, (t)"
label var cost_comply_cons "Fuel Cost, (USD)"
label var td_comply_cons "Damage AP2, (USD)"
label var td_inm_comply_cons "Damage, (USD)"
label var avg_speed "Speed, (km/hr)"
label var tkm_cons "Fuel, (t/km)"
label var tdtonne_pre_cons "Damage, (USD per t)"
label var tdtonne_inm_pre_cons "Damage, (USD per t)"

label var td_eca09_cons "Damage AP2, Comply Pre, (USD)"
label var td_inm_eca09_cons "Damage, Comply Pre, (USD)"
label var cost_eca09_cons "Fuel Cost, Comply Pre, (USD)"
label var td_eca11_cons "Damage AP2, Comply Pre, (USD)"
label var td_inm_eca11_cons "Damage, Comply Pre, (USD)"
label var cost_eca11_cons "Fuel Cost, Comply Pre, (USD)"

label var td_nox_cons "NOx Damages, (USD)"
label var td_nox_power "NOx Damage, (USD)"

label var dist_eca2009 "ECA Dist, (km)"
label var f_eca09_cons "ECA Fuel, (t)"
label var tkm_eca09_cons "ECA Fuel, (t/km)"
label var tkm_out09_cons "Out Fuel, (t/km)"


label var dist_eca2011 "ECA Dist, (km)"
label var f_eca11_cons "ECA Fuel, (t)"
label var tkm_eca11_cons "ECA Fuel, (t/km)"
label var tkm_out11_cons "Out Fuel, (t/km)"

label var hours "Travel Time, (hrs)"

label var built "Built (y)"
label var loa "Length (m)"
label var dwt "DWT (t)"
label var beam "Beam (m)"
label var draft "Draft"
label var power_kw "Power (kw)"

********************************************************************************

**** Classifying tracks as internal or as unclassified
gen same_port = ( portSN1==portSN2 ) 
gen single_port = ( portSN1==. | portSN2==.)
gen both_right = (right1==1 & right2==1) 

* number of right intersections
gen rights = 1
replace rights = 0 if right1~=1 & right2~=1
replace rights = 2 if right1==1 & right2==1 

gen internal = ( same_port==1 & both_right==1 ) | ( single_port==1 & rights==1 ) 
	* either enters port then exits, or enters and stays or exits and doesnt return
gen returner = ( same_port==1 & both_right==0 )
	* leaves then returns
gen unclassified = ( internal==0 ) & ( returner==0 )
	* all others
egen track_type = group(internal returner unclassified) , lab

*gen chk = internal + returner + unclassified		

* time series counts 
local sample eca_ind<=2 & main_sample==1 // & vesseltype_regstr=="Container"
lgraph dist my port_agg_i if `sample' & internal==1 & exposed==1, stat(count)

lgraph dist my port_agg_i if `sample' & unclassified==1 & exposed==1, stat(count)

* ignoring tracks that leave and return - small number
lgraph dist my port_agg_i if `sample' & returner==1 & exposed==1, stat(count)

		
**** GET # of tracks to expand counts in welfare by tracks that we can't complete
gen tmp_samp = (bad_months==0) & (speed_outlier==0) & interp_flag==0 & WCport==1 & vesseltype_regstr_filled~=""

preserve
collapse (count) dist if unclassified==1 & eca_ind==1 & tmp_samp==1 & port_agg_i~=. ///
	, by(port_agg vesseltype_regstr_filled)
rename dist tracks
outsheet using "rd_regs/track_counts_noroute.csv" , comma	replace	
restore

**** USE THIS TO GET # of internal tracks
preserve
collapse (count) dist if internal==1 & eca_ind==1 & tmp_samp==1 & exposed==1 ///
	, by(port_agg vesseltype_regstr_filled)
rename dist tracks
outsheet using "rd_regs/track_counts_internal.csv" , comma	replace	
restore

capture drop tmp_samp


/************************* RD for internal tracks******************/
levelsof eca_ind if eca_ind>0 , local(eca_inds)
foreach i of local eca_inds {
	qui sum date1 if eca_ind==`i'
	scalar eca`i'd = r(min)
}

*normalized time variable for RD
capture drop t_eca09 t_eca11
gen t_eca09 = date1 - eca1d
gen t_eca11 = date1 - eca2d


local run t_eca09
local donut 0
local controls built dwt loa draft beam power_kw USflag IFO LS
local xvar_rte "i.cut#c.`run'#i.route_id"
local bw 150
local vce "cluster ves_id"
local fe "i.route_id" // Using route name for this because it will uniquely identify ports for this sample
local folder "rd_regs"

capture drop cut
capture drop reg_samp
gen cut = eca1
label var cut "CA (1.5%), 2009" 
gen reg_samp = ( exposed==1 & main_sample==1 & internal==1 ) 

* looking at counts
*tab port_agg vesseltype_regstr if internal & main_sample & t_eca09<=150 & t_eca09>=-150

local report_extra_mean 0 

* getting mean fuel prices for predictions
* days pre ECA
local t_opts 30 0
foreach t of local t_opts {
	sum LS if `run'==-`t'
	local LSmean`t' = r(mean)
	sum IFO if `run'==-`t'
	local IFOmean`t' = r(mean)
}	

local yvars dist f_cons faux cost_comply_cons td_inm_comply_cons td_comply_cons // hours

levelsof port_agg if port_agg~="OtherWC", local(ports)
levelsof vesseltype_regstr, local(vestypes)
foreach prt of local ports {
	foreach vestype of local vestypes {
		local vestype2 : subinstr local vestype " " ""
		local prt1 : subinstr local prt " " ""
		local prt2 : subinstr local prt1 "/" ""
		
		local apprep replace
		
		foreach yvar of local yvars {

			local out_file "`folder'/rd_internal_`prt2'_2009_`vestype2'.tex"

			reghdfe `yvar' c.cut `xvar_rte' `controls' if reg_samp==1 & vesseltype_regstr=="`vestype'" & port_agg=="`prt'" /// 
					& ( `run'<=`bw' & `run'>=-`bw' ) & ( `run'<-`donut' | `run'>`donut' ) , absorb(`fe') vce(`vce')
			estimates store regresults
			
			* predicting outcome at t=0 and t=-30
			margins if e(sample) & cut==0 , predict(xb) at( `run'=( -30 ) LS=(`LSmean30') IFO=(`IFOmean30') ) /// 
														  at( `run'=( 0 ) LS=(`LSmean0') IFO=(`IFOmean0') )  post

			local mean_t30 = _b[1bn._at]
			local mean_t0 = _b[2._at]	
						
			if `report_extra_mean'==1 {
				local mean_text ", Mean (t=0) , `mean_t0', Mean (t=-30), `mean_t30'"
			}
			else {
				local mean_text ", Mean (t=0) , `mean_t0'"
			}
						
			estimates restore regresults
			
			*sum `yvar' if e(sample) & cut==0 
			*local mean = r(mean)
			
			outreg2 using `out_file', `apprep' lab keep( cut ) tex(frag) nonotes nocon /// 
			addstat(Vessels , e(N_clust) `mean_text' )  
			local apprep append
			
		}
	}
}

